{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4da87a7d",
   "metadata": {},
   "source": [
    "# Python机器学习（第15期）第2课书面作业\n",
    "1. 用韦来生书第83页定义的后验众数估计，后验中位数估计，后验期望估计等方法，在原样例代码基础之上，写一小段python代码计算短信计数的例子中lambda1，lambda2，tau等几项参数的点估计\n",
    "2. 韦来生书第75页第2题（可选）\n",
    "3. 通过计算证明伽马分布和泊松分布是共轭分布（参考韦来生书第73页内容）\n",
    "## 第1题\n",
    "韦来生书中应该是在81页（不是83页）定义了贝叶斯点估计方法，如后验众数估计，后验中位数估计，后验期望估计。  \n",
    "定义 4.1.1 用后验密度 $\\pi(\\theta|x)$ 达到最大值时$\\theta$之值作为估计量，称为$\\theta$的后验众数估计或广义最大似然估计，记为$\\hat{\\theta}_{MD}$。用后验分布的中位数作为$\\theta$的估计量，称为后验中位数估计，记为$\\hat{\\theta}_{ME}$。用后验分布的期望值作为$\\theta$的估计量，称为后验期望估计，记为$\\hat{\\theta}_{E}$。  \n",
    "下面针对《Probabilistic Programming and Bayesian Methods for Hackers》书中的计算短信计数的例子，计算其中lambda1，lambda2，tau等几项参数的点估计。  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a24708ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAukAAAD7CAYAAADXRQqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmzklEQVR4nO3deZxkdXnv8c9XVkEQkIEMmwRFDHoBdSQqXkNYcsEFuEZQImZEFDUkkXsVBOO9USMKmhjXGHHBEVAvEhE0uEwGwbgzGFYHMkhkkWEYUXYDDDz3j3M6FG339Olmqru66/N+vepVdbbfeepXVd1P/eo556SqkCRJkjQ4HjPTAUiSJEl6JJN0SZIkacCYpEuSJEkDxiRdkiRJGjAm6ZIkSdKAMUmXJEmSBoxJujTEkvxjkv+zhuWV5MmTaG9S62tsSd6W5FN9anvKr1GSnyfZb5xl/z3JNY8uugn3/+ok3+3nPuaaJK9M8q2ZjkPS5JmkS3NUm1D9JsldSW5P8v0kb0jyX5/7qnpDVf3NFNu/MMlr117EM6N9HntPcpu9k9y0lvb/W21V1Xuqalb1bVX9a1XtMjK9poRe/ZFkx/ZL2Loj86rqzKr6oz7sa7Mkn0tyc5JfJflIkqzt/UjDzCRdmtteUlWbAE8ETgbeCnx6ZkMaDGn4N1CzUm8iPkO2AP4N2LW9vRg4dEYjkuYY/0FJQ6Cq7qiq84CXAwuTPB0gyWeTvHtkvSTHJVnRjo69Zrz2kpwE/Hfgo0nuTvLRnsX7JVme5NdJPtY7upbkNUmWtcu+meSJ47T/W6PLvSOzSfZMsjTJnUlWJvlAz3rPaX81uD3JZb2j5O2o+UlJvgfcC+w0ah/jttuzzsbA14Ft2ud+d5JtkjwmyQlJfpbktiRnJdmi3ebjSc7uaeOUJEvW0NY7kpzRrjsyOrowyQ1Jfpnkr3raemySRW2fLktyfIdR/jFfoyRPSnJBG/8vk5yZZLNR2z47yU/bbU9LsuHo1yzJ6cAOwFfb53R8kg2TnNG2fXuSi5NsPVZwSbZP8uUkq9r1Pzpq+d+2+/+PJAf2zD+y7YO7klyX5PU9y/ZOclOSNye5tX2fH9mz/AlJvtq+9hcneXd6SmuSPDXJ4jSjxtckOWy8zm1fw/Pada9N8rqe+b8ZeV+0857R9vV67fS4n5H2fXBMkuXA8jF2/Z32/va235+bUSVCbRt/1r7+dyX5m/Z1/0H73M9Ksn7P+i9Ocmke/jVuN4Cquq6q/r6qbq+qW4B/B8Z8PSVNUVV58+ZtDt6AnwP7jTH/BuCN7ePPAu9uHx8ArASeDmwMfB4o4MnjtH8h8NpR8wr4GrAZTZK2CjigXXYIcC3we8C6wNuB74/T9t7ATeM9H+AHwKvax48DntM+3ha4DXghzSDE/u30vJ6YbwCe1saw3qh9jNlux/iOBX4IbAdsAHwC+EK7bCOaJObVNF9ufglst4a23gGc0T7ese3XTwKPBXYH7gN+r11+MnARsHm778tHtzeJ1+jJbZ9tAMyjSfo+OOo1uBLYnmYk9Xs8/P55xPNg1PsPeD3w1bYv1gGeBWw6RnzrAJcBf0/zPtwQeH677NXAA8Dr2vXeCNwMpF3+IuBJQIA/oPki9sye+FYD7wLWo3mP3Ats3i7/YnvbiGZk+Ebgu+2yjdvpI2neN89sX8OnjdPHFwH/0Ma+R9vH+7bLLgBe17Pu+4F/7PIZaV+7xW3fP3aM/Y68V9btmffqkefR08Z5wKY0n4P7gCU0X1gfD/wUWNiu+0zgVuD32/5e2L6uG4za78uAXwHbz/TfPW/e5tLNkXRp+NxM809+tMOA06rqyqq6hyZRnIqTqxlduwH4Nk2SAk2S9t6qWlZVq4H3AHtknNH0CTwAPDnJllV1d1X9sJ1/BHB+VZ1fVQ9V1WJgKU1CNuKzVXVVVa2uqgc6ttvF64G/qqqbquo+mv57WZJ1q+reNrYPAGcAf1FVk61pf2dV/aaqLqNJYndv5x8GvKeqft22+eEObY35GlXVtVW1uKruq6pVbbx/MGrbj1bVjVX1K+Ak4PCO8T8APIHmS9+DVXVJVd05xnp7AtsAx1XVPVX1n1XVe7Do9VX1yap6EFgEzKcdwa2qf66qn1XjIuBbNF+KemN4V1U9UFXnA3cDuyRZB/hj4K+r6t6q+mnb9ogXAz+vqtPa981PgH+iSU4fIcn2wPOBt7axXwp8CnhVu8rnR/qs/QXjFe086PYZeW9V/aqqfjNWJ3d0SlXdWVVX0Xzp+lY1I+N30Pyy84x2vdcBn6iqH7Wv2SKapP45Pc/3+TRfIA+qqhsfRUySRjFJl4bPtjSjXqNtQzNaOOL6KbZ/S8/je2lGpKGpi/9Q+7P57W0MaeOZrKOApwBXt6UJL+7Zx6Ej+2j383yaRG7EmhKJ8drt4onAOT37XQY8yMMJ5I+B62ie81mTaHfEeP06+nXrkiiN2VaSrZJ8MckvktxJ84Viy1Hbjn6PbNNhfwCnA98EvpimnOp9IyUeo2xPk4ivnij29ssPPfEfmOSHbZnJ7TRfznrjv21UuyPPfR7NyPV4/fhE4PdHva9eCfzOGPFtA/yqqu7qmXc9D7/Pzwaem2Qb4AU0I9v/2rOfiT4jayMRXtnz+DdjTPd+Zt886nlvzyNf8z8D/n7UFylJa4FJujREkjyb5h/+WP9QV9D8Ax6xwwTN1SR3fyPw+qrarOf22Kr6/hjr3kNTdjAS9zo0iVSz46rlVXU4sBVwCnB2mvruG4HTR+1j46o6uUvca2j3t1Yd5/kdOGrfG1bVL9rncAxNGcnNwPFd4uloBU2Zy4jtx1uxg/fSxLNbVW1KM/o/+owdo98jN4/T1iOeVzt6/c6q2hV4Hs3o9J+Osd2NwA6Z5IGRSTagGd3+W2DrqtoMOH+M+MeyiqYUZrx+vBG4aNRr+7iqeuMYbd0MbJFkk555OwC/AKiq22lG+A8D/oSmJGqkr7p8Rtb0fnm076XRbgROGhXPRlX1hZ515jP+e0DSo2CSLg2BJJu2o8JfpKl1vmKM1c4CXp1k1yQbAX89QbMrGXXg5QT+ETgxydPamB6fZLyzQfw7sGGSF7WjrW+nSXBHns8RSeZV1UPA7e3sB2lGfl+S5H8kWSfNwYp7J9lu9A7GsoZ2R1sJPCHJ40c9v5NGShOSzEtycPv4KcC7aZLeVwHHJ9ljDW1Nxlk0/bp5km2BP59iOwCb0JSA3N62ddwY6xyTZLv24Me3Af9vnLYe8f5I8odJ/lv7hetOmtKTsfr2xzRfPE5OsnH7Gu7VIfb1ad4jq4DVaQ4o7XTqwbZ05svAO5JslOSpPPILxNeApyR5VZL12tuzk/zeGG3dCHwfeG8b+240v9Cc2bPa59v2/5iHS11gcp+RsawCHmJyn8s1+STwhiS/n8bG7Wey9wvIH/PI5yZpLTFJl+a2rya5i2ZE7K9oaoyPHGvFqvo68EGaA9uube/X5EM0Nde/TjJhHXRVnUMzOv3FtpTiSuDAcda9g+Zn9E/RjEDeA/TWcB8AXJXk7jaOV7T1vzcCB9Mkj6tonvdxdP9bN2a7Y8R3NfAF4Lq2DGCbdv3zgG+1ff5DmhKJdWm+PJxSVZdV1fI2vtOTbDBOW5PxLpq++Q/gX2jKKe6bZBsj3klzsOAdwD/TJK6jfZ5mJPi69vbuMdaBZlT+7e1zegtNacjZNAn6MpqDK88YvVGbML+E5iDWG9rn9vKJAm/LS/6S5kvLr2lGqc+baLsef05z4OQtNKU5X6Dtx7btP6KpH7+5XecUer44jnI4zUGcNwPn0NS6L+5Zfh6wM7CyPcZg5Dl0/oyMpS3/OQn4Xtvvz5lomwnaW0pTl/5Rmj69luZA1F5n0vwqIGktGzkiXpI0ByR5I82Xi9EHfGoSkpwC/E5VLZzpWCQNJ0fSJWkWSzI/yV5pztO+C/BmmtFbTUKa86Dv1pZ17ElTomI/SpoxM33FMknSo7M+zTnZf5emjv6LNOfo1uRsQlPisg3NucH/Djh3RiOSNNQsd5EkSZIGjOUukiRJ0oAxSZckSZIGzKyoSd9yyy1rxx13nOkwJEmSpLXmkksu+WVVzRtr2axI0nfccUeWLl0602FIkiRJa02S68dbZrmLJEmSNGBM0iVJkqQBY5IuSZIkDRiTdEmSJGnAmKRLkiRJA8YkXZIkSRowJumSJEnSgJkV50mXtHY9+/0XrHH5xcftM02RSJKksfRtJD3JLkku7bndmeTYJFskWZxkeXu/eb9ikCRJkmajviXpVXVNVe1RVXsAzwLuBc4BTgCWVNXOwJJ2WpIkSVJrumrS9wV+VlXXAwcDi9r5i4BDpikGSZIkaVaYriT9FcAX2sdbV9UKgPZ+q2mKQZIkSZoV+n7gaJL1gYOAEye53dHA0QA77LBDHyKTJGlyPOha0nSZjpH0A4GfVNXKdnplkvkA7f2tY21UVadW1YKqWjBv3rxpCFOSJEkaDNORpB/Ow6UuAOcBC9vHC4FzpyEGSZIkadboa5KeZCNgf+DLPbNPBvZPsrxddnI/Y5AkSZJmm77WpFfVvcATRs27jeZsL5IkSZLGMF1nd5EkSZLUkUm6JEmSNGBM0iVJkqQBY5IuSZIkDRiTdEmSJGnAmKRLkiRJA8YkXZIkSRowJumSJEnSgBn3YkZJrgBqvOVVtVtfIpIkSZKG3JquOPri9v6Y9v709v6VwL19i0iSJEkacuMm6VV1PUCSvapqr55FJyT5HvCufgcnSZIkDaMuNekbJ3n+yESS5wEb9y8kSZIkabitqdxlxFHAZ5I8nqZG/Q7gNX2NSpIkSRpiEybpVXUJsHuSTYFU1R39D0uSJEkaXhOWuyTZOsmngf9XVXck2TXJUdMQmyRJkjSUutSkfxb4JrBNO/3vwLF9ikeSJEkael2S9C2r6izgIYCqWg082KXxJJslOTvJ1UmWJXluki2SLE6yvL3f/FHEL0mSJM05XZL0e5I8gfbCRkmeQ3PwaBcfAr5RVU8FdgeWAScAS6pqZ2BJOy1JkiSp1eXsLm8GzgOe1J4ffR7wsok2ag80fQHwaoCquh+4P8nBwN7taouAC4G3TjJuSZIkac7qdHaXJH8A7AIEuKaqHujQ9k7AKuC0JLsDlwBvArauqhVt2yuSbDXl6CVJkqQ5qMvZXS4Djgf+s6qu7JigQ/MF4JnAx6vqGcA9TKK0JcnRSZYmWbpq1aqum0mSJEmzXpea9IOA1cBZSS5O8pYkO3TY7ibgpqr6UTt9Nk3SvjLJfID2/taxNq6qU6tqQVUtmDdvXofdSZIkSXPDhEl6VV1fVe+rqmcBfwLsBvxHh+1uAW5Msks7a1/gpzT17QvbeQuBc6cSuCRJkjRXdTlwlCQ7AocBL6c5/eLxHdv/C+DMJOsD1wFH0nwxOKu9ININwKGTjFmSJEma0yZM0pP8CFgP+BJwaFVd17XxqroUWDDGon27tiFJkiQNmy4j6Qur6uq+RyJJkiQJWEOSnuSIqjoDeGGSF45eXlUf6GtkkiRJ0pBa00j6xu39JtMRiIbTs99/wbjLLj5un2mMRJIkaXCMm6RX1Sfa+3dOXziSJEmSulzM6ClJliS5sp3eLcnb+x+aJEmSNJy6XMzok8CJwAMAVXU58Ip+BiVJkiQNsy5J+kZV9eNR81b3IxhJkiRJ3ZL0XyZ5ElAASV4GrOhrVJIkSdIQ63Ke9GOAU4GnJvkF8B/AEX2NSpIkSRpiEybp7RVG90uyMfCYqrqr/2FJkiRJw6vL2V3ek2Szqrqnqu5KsnmSd09HcJIkSdIw6lKTfmBV3T4yUVW/Bn7rCqSSJEmS1o4uNenrJNmgqu4DSPJYYIP+hiVJksArM0vDqkuSfgawJMlpNGd4eQ2wqK9RSZIkSUOsy4Gj70tyObAfEOBvquqbfY9MkiRJGlJdRtIBlgGrq+pfkmyUZBPP8iJJkiT1R5ezu7wOOBv4RDtrW+ArfYxJkiRJGmpdzu5yDLAXcCdAVS0HturSeJKfJ7kiyaVJlrbztkiyOMny9n7zqQYvSZIkzUVdkvT7qur+kYkk69IcQNrVH1bVHlW1oJ0+AVhSVTsDS9ppSZIkSa0uSfpFSd4GPDbJ/sCXgK8+in0ezMNnh1kEHPIo2pIkSZLmnC5J+luBVcAVwOuB84G3d2y/gG8luSTJ0e28ratqBUB7P2bpTJKjkyxNsnTVqlUddydJkiTNfms8u0uSxwCXV9XTgU9Oof29qurmJFsBi5Nc3XXDqjoVOBVgwYIFkymvkSRJkma1NSbpVfVQksuS7FBVN0y28aq6ub2/Nck5wJ7AyiTzq2pFkvnArVOKXJKkAeVVQiU9Wl3KXeYDVyVZkuS8kdtEGyXZOMkmI4+BPwKuBM4DFrarLQTOnVrokiRJ0tzU5WJG75xi21sD5yQZ2c/nq+obSS4GzkpyFHADcOgU25ckSZLmpAmT9Kq6aCoNV9V1wO5jzL8N2HcqbUqSJEnDoEu5iyRJkqRpZJIuSZIkDZhJJelJNk+yW7+CkSRJktQhSU9yYZJNk2wBXAacluQD/Q9NkiRJGk5dRtIfX1V3Ai8FTquqZwH79TcsSZIkaXh1SdLXbS86dBjwtT7HI0mSJA29Lkn6u4BvAtdW1cVJdgKW9zcsSZIkaXh1OU/6l4Av9UxfB/xxP4OSJEmShtmESXqSD48x+w5gaVWdu/ZDkiRJkoZbl3KXDYE9aEpclgO7AVsARyX5YN8ikyRJkobUhCPpwJOBfapqNUCSjwPfAvYHruhjbJIkSdJQ6jKSvi2wcc/0xsA2VfUgcF9fopIkSZKGWJeR9PcBlya5EAjwAuA9STYG/qWPsUmSJElDqcvZXT6d5HxgT5ok/W1VdXO7+Lh+BidJkiQNownLXZLsBdzZnsnlccBbkjyx75FJkiRJQ6pLTfrHgXuT7E4zcn498Lm+RiVJkiQNsS5J+uqqKuBg4MNV9SFgk647SLJOkn9L8rV2eoski5Msb+83n1rokiRJ0tzUJUm/K8mJwBHAPydZB1hvEvt4E7CsZ/oEYElV7QwsaaclSZIktbok6S+nOdXiUVV1C80pGd/fpfEk2wEvAj7VM/tgYFH7eBFwSNdgJUmSpGEwYZJeVbdU1QeAy5JsAdwNfK1j+x8Ejgce6pm3dVWtaNteAWw11oZJjk6yNMnSVatWddydJEmSNPt1ObvL65OsBC4HLmlvSzts92Lg1qq6ZCqBVdWpVbWgqhbMmzdvKk1IkiRJs1KXixm9BXhaVf1ykm3vBRyU5IXAhsCmSc4AViaZX1UrkswHbp1ku5IkSdKc1iVJ/xlw72QbrqoTgRMBkuwNvKWqjkjyfmAhcHJ7f+5k25YkSdLc8+z3X7DG5Rcft880RTLzuiTpJwLfT/IjmgNIAaiqv5ziPk8GzkpyFHADcOgU25EkSZLmpC5J+ieAC4AreOQBoJ1V1YXAhe3j24B9p9KOJEmSNAy6JOmrq+p/9z0SSZIkSUC386R/uz0d4vz2aqFbtKdilCRJktQHXUbS/6S9P7FnXgE7rf1wJEmSpDVb0wGmc+Xg0gmT9Kr63ekIRJIkSVKjS7nLf0lyar8CkSRJktSYVJIOLOhLFJIkSZL+y2STdK8OKkmSJPXZhDXpSQ6tqi8BVNUBo+fpkbxSliRJg8v/05otuoykn9hxniRJkqS1YNyR9CQHAi8Etk3y4Z5FmwKr+x2YJEmSNKzWVO5yM7AUOAi4pGf+XcD/6mdQkiRJ0jAbN0mvqsuAy5Isqaqbepcl2QX4db+DkyRJkoZRl5r0JUkOG5lI8mbgnP6FJEmSJA23Cc/uAuwNnJrkUGBrYBmwZz+DkiRJkobZhEl6Va1I8g2aM7o8BJxYVXf3PTJJM8rTlEmSNHO6nCd9MbACeDqwHfCZJN+pqrf0OzhJkiRpGHWpSf9YVf1pVd1eVVcCzwXumGijJBsm+XGSy5JcleSd7fwtkixOsry93/xRPgdJkiRpTulS7vKVJM8Hdq6q04DNgTM6tH0fsE9V3Z1kPeC7Sb4OvBRYUlUnJzkBOAF469SfgiRJWlvWVOpmmdvsYLni3DDhSHqSv6ZJokeuMro+HZL0aozUrq/X3go4GFjUzl8EHDK5kCVJkqS5rcvZXf4n8AzgJwBVdXOSTbo0nmQdmgshPZmmbOZHSbauqhVtWyuSbDXOtkcDRwPssMMOXXYnSZKkDhxtH3xdatLvr6qiGQUnycZdG6+qB6tqD5oDTvdM8vRJbHtqVS2oqgXz5s3rupkkSZI063VJ0s9K8glgsySvA/4F+NRkdlJVtwMXAgcAK5PMB2jvb51MW5IkSdJc1+XA0b9Nsj9wJ7AL8H+ravFE2yWZBzxQVbcneSywH3AKcB6wEDi5vT/3UcQvTYo/70mSpNmgy3nST6mqtwKLx5i3JvOBRW1d+mOAs6rqa0l+QDM6fxRwA3Do1MOXJEmS5p4uB47uz2+fIvHAMeY9QlVdTnPA6ej5twH7dg1QkiRJGjbjJulJ3gj8GbBTkst7Fm0CfK/fgUmSJEnDak0j6Z8Hvg68l+aCQyPuqqpf9TUqSZIkaYiNm6RX1R3AHcDh0xeOJE0fr6woSY+eJ2Xojy6nYJQkSZI0jcZN0pNsMJ2BSJIkSWqsaST9BwBJTp+mWCRJkiSx5gNH10+yEHhekpeOXlhVX+5fWJIkSdLwWlOS/gbglcBmwEtGLSvAJF2SJEnqgzWd3eW7wHeTLK2qT09jTJIkSdJQ63LF0dOT/CXwgnb6IuAfq+qB/oXVzdo6fZqnDpIkSdIg6ZKk/wOwXnsP8Crg48Br+xWUJEmSNMy6JOnPrqrde6YvSHJZvwKSJEmShl2XJP3BJE+qqp8BJNkJeLC/YUnSo+PVRCX1m+Wy6qcuSfpxwLeTXAcEeCJwZF+jkiRJkobYhEl6VS1JsjOwC02SfnVV3df3yDTtuo4IOEK59tmnkh4NR3Slmbe2P4ddRtJpk/LLJ9WyJEmSpCl5TL8aTrJ9km8nWZbkqiRvaudvkWRxkuXt/eb9ikGSJEmajTqNpE/RauDNVfWTJJsAlyRZDLwaWFJVJyc5ATgBeGsf45AkaehZEiPNLhOOpCdZ0mXeaFW1oqp+0j6+C1gGbAscDCxqV1sEHDKJeCVJkqQ5b9yR9CQbAhsBW7YlKWkXbQpsM5mdJNkReAbwI2DrqloBTSKfZKtxtjkaOBpghx12mMzuJGmtcfRRGgweYK9hs6Zyl9cDx9Ik5JfwcJJ+J/CxrjtI8jjgn4Bjq+rOJBNtAkBVnQqcCrBgwYLquj9JkiRpths3Sa+qDwEfSvIXVfWRqTSeZD2aBP3MqvpyO3tlkvntKPp84NaptC1JkiTNVV3Ok/6RJM8Dduxdv6o+t6bt0gyZfxpYVlUf6Fl0HrAQOLm9P3fyYUu/bbp/CvWn15lh+Yn6xc90d34Opf6bMElPcjrwJOBS4MF2dgFrTNKBvYBXAVckubSd9zaa5PysJEcBNwCHTjpqSZIkaQ7rcgrGBcCuVTWpuvCq+i4P17GPtu9k2pIkSZKGSZck/Urgd4AVfY5F0lrgT/aSZhNLZ6SxdUnStwR+muTHwH0jM6vqoL5FJUmSJA2xLkn6O/odhCQNA0cMpdljbf4q6S+cmoouZ3e5aDoCkSRJktTocnaXu2jO5gKwPrAecE9VbdrPwCRJkqRh1WUkfZPe6SSHAHv2KyBJ0uCxVEeSptdjJrtBVX0F8K+xJEmS1Cddyl1e2jP5GJrzpk/qnOnqLw9IkTQoBunKvyP79FcASbNRl7O7vKTn8Wrg58DBfYlGkiRJUqea9COnIxBJkiRJjS7lLtsBHwH2oilz+S7wpqq6qc+xzWmWqEiaDv6tkaTxDXI5XJcDR08DzgO2AbYFvtrOkyRJktQHXWrS51VVb1L+2STH9imeMS275a5HfNNx9Eeamxz1bdgPkgaFf48aM9EPXUbSf5nkiCTrtLcjgNv6Eo0kSZKkTkn6a4DDgFuAFcDL2nmSJEmS+qDL2V1uAA6ahlikOWdtHpAyyAe3SNIg8O+k5pIJR9KTLEqyWc/05kk+02G7zyS5NcmVPfO2SLI4yfL2fvMpRy5JkiTNUV0OHN2tqm4fmaiqXyd5RoftPgt8FPhcz7wTgCVVdXKSE9rpt3YPd7jM5lHY2Ry7JGnt8KBDaeq61KQ/pnfEO8kWdCuT+Q7wq1GzDwYWtY8XAYd0C1OSJEkaHl1G0v8O+H6Ss2kuZnQYcNIU97d1Va0AqKoVSbaaYjuSJEnSnNVlRPxzSZYC+wABXlpVP+13YEmOBo4GWH8zc3lJ6mUZmCTNbV1G0mmT8rWRmK9MMr8dRZ8P3LqGfZ4KnAqw8Xa71FrYtyRJkjQrdErS16LzgIXAye39udO8f0lrmQeGSZK09nU5cHRKknwB+AGwS5KbkhxFk5zvn2Q5sH87LUmSJKlH30bSq+rwcRbt2699SpIkSXPBdJe7TDsPrpIkSdJs07dyF0mSJElTM+dH0tcmD5CTpsZftNTL90PDfpC0Jo6kS5IkSQPGJF2SJEkaMJa7SKP4E7SmYjaXw83m2CVprnIkXZIkSRowjqRLkiRpzpntv4w7ki5JkiQNGJN0SZIkacBY7iJJkjTDZntphtY+R9IlSZKkAWOSLkmSJA0Yk3RJkiRpwJikS5IkSQPGA0clSQPJK6FKv83PxfCYkZH0JAckuSbJtUlOmIkYJEmSpEE17Ul6knWAjwEHArsChyfZdbrjkCRJkgbVTIyk7wlcW1XXVdX9wBeBg2cgDkmSJGkgzUSSvi1wY8/0Te08SZIkSUCqanp3mBwK/I+qem07/Spgz6r6i1HrHQ0c3U4+HbhyWgMVwJbAL2c6iCFkv88c+35m2O8zw36fGfb7zBnEvn9iVc0ba8FMnN3lJmD7nuntgJtHr1RVpwKnAiRZWlULpic8jbDfZ4b9PnPs+5lhv88M+31m2O8zZ7b1/UyUu1wM7Jzkd5OsD7wCOG8G4pAkSZIG0rSPpFfV6iR/DnwTWAf4TFVdNd1xSJIkSYNqRi5mVFXnA+dPYpNT+xWL1sh+nxn2+8yx72eG/T4z7PeZYb/PnFnV99N+4KgkSZKkNZuRK45KkiRJGt9AJ+lJDkhyTZJrk5ww0/HMZUk+k+TWJFf2zNsiyeIky9v7zWcyxrkoyfZJvp1kWZKrkrypnW/f91GSDZP8OMllbb+/s51vv0+DJOsk+bckX2un7fdpkOTnSa5IcmmSpe08+77PkmyW5OwkV7d/659rv/dXkl3a9/nI7c4kx862fh/YJD3JOsDHgAOBXYHDk+w6s1HNaZ8FDhg17wRgSVXtDCxpp7V2rQbeXFW/BzwHOKZ9n9v3/XUfsE9V7Q7sARyQ5DnY79PlTcCynmn7ffr8YVXt0XMaOvu+/z4EfKOqngrsTvPet9/7qKquad/newDPAu4FzmGW9fvAJunAnsC1VXVdVd0PfBE4eIZjmrOq6jvAr0bNPhhY1D5eBBwynTENg6paUVU/aR/fRfPHe1vs+76qxt3t5HrtrbDf+y7JdsCLgE/1zLbfZ45930dJNgVeAHwaoKrur6rbsd+n077Az6rqemZZvw9ykr4tcGPP9E3tPE2fratqBTTJJLDVDMczpyXZEXgG8CPs+75rSy4uBW4FFleV/T49PggcDzzUM89+nx4FfCvJJe1VvcG+77edgFXAaW2J16eSbIz9Pp1eAXyhfTyr+n2Qk/SMMc9T0WhOSvI44J+AY6vqzpmOZxhU1YPtT6HbAXsmefoMhzTnJXkxcGtVXTLTsQypvarqmTRlpMckecFMBzQE1gWeCXy8qp4B3MOAl1jMJe1FMw8CvjTTsUzFICfpNwHb90xvB9w8Q7EMq5VJ5gO097fOcDxzUpL1aBL0M6vqy+1s+36atD89X0hzTIb93l97AQcl+TlNCeM+Sc7Afp8WVXVze38rTX3untj3/XYTcFP7Sx3A2TRJu/0+PQ4EflJVK9vpWdXvg5ykXwzsnOR3229CrwDOm+GYhs15wML28ULg3BmMZU5KEppaxWVV9YGeRfZ9HyWZl2Sz9vFjgf2Aq7Hf+6qqTqyq7apqR5q/6RdU1RHY732XZOMkm4w8Bv4IuBL7vq+q6hbgxiS7tLP2BX6K/T5dDufhUheYZf0+0BczSvJCmvrFdYDPVNVJMxvR3JXkC8DewJbASuCvga8AZwE7ADcAh1bV6INL9SgkeT7wr8AVPFyj+zaaunT7vk+S7EZz0NA6NIMVZ1XVu5I8Aft9WiTZG3hLVb3Yfu+/JDvRjJ5DU4Lx+ao6yb7vvyR70BwovT5wHXAk7d8d7Pe+SbIRzbGNO1XVHe28WfV+H+gkXZIkSRpGg1zuIkmSJA0lk3RJkiRpwJikS5IkSQPGJF2SJEkaMCbpkiRJ0oAxSZekAZXkCUkubW+3JPlF+/juJP/Qp30em+RPx5i/Y5Ir1+J+1k/ynSTrrq02JWku8Y+jJA2oqroN2AMgyTuAu6vqb/u1vzZhfg3NFRH7qqruT7IEeDlwZr/3J0mzjSPpkjTLJNk7ydfax+9IsijJt5L8PMlLk7wvyRVJvpFkvXa9ZyW5KMklSb45cmnsUfahuYT26p5tLkvyA+CYnv3vmORfk/ykvT2vnX96koN71jszyUFJnpbkx+2vAJcn2bld5SvAK/vRR5I025mkS9Ls9yTgRcDBwBnAt6vqvwG/AV7UJuofAV5WVc8CPgOMdQXnvYBLeqZPA/6yqp47ar1bgf2r6pk0I+Efbud/iuZqiiR5PPA84HzgDcCHqmoPYAFwU7v+lcCzp/icJWlOs9xFkma/r1fVA0muANYBvtHOvwLYEdgFeDqwOAntOivGaGc+sAz+K8nerKouapedDhzYPl4P+Gh7ufMHgacAVNVFST6WZCvgpcA/VdXqdiT+r5JsB3y5qpa36z+Y5P4km1TVXWupLyRpTjBJl6TZ7z6AqnooyQNVVe38h2j+zge4aowR8dF+A2zYPg5Q46z3v4CVwO40v8j+Z8+y02lKWF5BU99OVX0+yY9oRvu/meS1VXVBu/4Go7aXJGG5iyQNg2uAeUmeC5BkvSRPG2O9ZcCTAarqduCOJM9vl/XWjj8eWFFVDwGvohmZH/FZ4Ni2java/e0EXFdVHwbOA3Zr5z8BWFVVDzz6pyhJc4tJuiTNcVV1P/Ay4JQklwGX0tSLj/Z14AU900cCH2vLVX7TM/8fgIVJfkhT6nJPz75W0iT7p/Ws/3LgyiSXAk8FPtfO/0OamnVJ0ih5+FdRSdKwS3IOcPxI3fgUtt+Iphb+mVV1xwTrfhk4saqumcq+JGkucyRdktTrBJoDSCctyX7A1cBHOiTo6wNfMUGXpLE5ki5JkiQNGEfSJUmSpAFjki5JkiQNGJN0SZIkacCYpEuSJEkDxiRdkiRJGjAm6ZIkSdKA+f8Sa18Bg6y/4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 900x252 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from IPython.core.pylabtools import figsize\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "figsize(12.5, 3.5)\n",
    "count_data = np.loadtxt(\"data/txtdata.csv\")\n",
    "n_count_data = len(count_data)\n",
    "plt.bar(np.arange(n_count_data), count_data, color=\"#348ABD\")\n",
    "plt.xlabel(\"Time (days)\")\n",
    "plt.ylabel(\"count of text-msgs received\")\n",
    "plt.title(\"Did the user's texting habits change over time?\")\n",
    "plt.xlim(0, n_count_data);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fa90d81b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "\n",
    "with pm.Model() as model:\n",
    "    alpha = 1.0/count_data.mean()  # Recall count_data is the\n",
    "                                   # variable that holds our txt counts\n",
    "    lambda_1 = pm.Exponential(\"lambda_1\", alpha)\n",
    "    lambda_2 = pm.Exponential(\"lambda_2\", alpha)\n",
    "    \n",
    "    tau = pm.DiscreteUniform(\"tau\", lower=0, upper=n_count_data - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ed5a8539",
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    idx = np.arange(n_count_data) # Index\n",
    "    lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f016c468",
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    observation = pm.Poisson(\"obs\", lambda_, observed=count_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0b3c6626",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-7-2bd7fdec376c>:4: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n",
      "  trace = pm.sample(10000, tune=5000,step=step)\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">Metropolis: [tau]\n",
      ">Metropolis: [lambda_2]\n",
      ">Metropolis: [lambda_1]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='30000' class='' max='30000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [30000/30000 13:57<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 5_000 tune and 10_000 draw iterations (10_000 + 20_000 draws total) took 850 seconds.\n",
      "The number of effective samples is smaller than 25% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "### Mysterious code to be explained in Chapter 3.\n",
    "with model:\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(10000, tune=5000,step=step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5fdc07cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_1_samples = trace['lambda_1']\n",
    "lambda_2_samples = trace['lambda_2']\n",
    "tau_samples = trace['tau']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33c4b8ed",
   "metadata": {},
   "source": [
    "下面分别计算lambda_1、lambda_2、tau的后验众数估计，后验中位数估计，后验期望估计："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "6ba6b177",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "1ccfcb8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                  Theta_MD Theta_ME Theta_E\n",
      "lambda_1_samples: 17.7936, 17.7505, 17.7478\n",
      "lambda_2_samples: 22.4244, 22.7039, 22.7273\n",
      "     tau_samples: 45.0000, 44.0000, 44.3024\n"
     ]
    }
   ],
   "source": [
    "print('                  Theta_MD Theta_ME Theta_E')\n",
    "print('lambda_1_samples: %.4f, %.4f, %.4f'%(stats.mode(lambda_1_samples)[0][0],\n",
    "                                            np.median(lambda_1_samples),\n",
    "                                            np.mean(lambda_1_samples)))\n",
    "print('lambda_2_samples: %.4f, %.4f, %.4f'%(stats.mode(lambda_2_samples)[0][0],\n",
    "                                            np.median(lambda_2_samples),\n",
    "                                            np.mean(lambda_2_samples)))\n",
    "print('     tau_samples: %.4f, %.4f, %.4f'%(stats.mode(tau_samples)[0][0],\n",
    "                                            np.median(tau_samples),\n",
    "                                            np.mean(tau_samples)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e363991",
   "metadata": {},
   "source": [
    "上面输出中Theta_MD,Theta_ME和Theta_E分别表示：后验众数估计，后验中位数估计，后验期望估计。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8aac3a7",
   "metadata": {},
   "source": [
    "## 第2题\n",
    "韦来生书第75页第2题：\n",
    "设某校学生的身高（单位：cm）服从$N(\\theta,5^2)$，今从该校学生中随机抽取10人测量其身高，其平均高度为175.34cm，设平均高度$\\theta$的先验分布为$N(172.72,2.56)$，求$\\theta$的后验分布。\n",
    "\n",
    "答：  \n",
    "从题目可以看出是属于”当$\\sigma^2$已知时，均值参数$\\theta$的后验分布“情况，按书中3.2.16及3.2.17式可见：\n",
    "$$\n",
    "\\pi(\\theta|\\overline{x})=N(\\mu_n(\\overline{x}),\\eta_n^2) \\\\\n",
    "$$\n",
    "其中上式的部分符号定义如下：\n",
    "$$\n",
    "\\mu_n(\\overline{x})=\\frac{\\sigma_n^2\\mu+\\tau^2\\overline{x}}{\\sigma_n^2+\\tau^2} \\\\\n",
    "\\eta_n^2=\\frac{\\sigma_n^2\\tau^2}{\\sigma_n^2+\\tau^2} \\\\\n",
    "\\sigma_n^2=\\sigma^2/n \\\\\n",
    "$$\n",
    "可知：\n",
    "\n",
    "1. $\\sigma^2=5^2, n=10$，因此$\\sigma_n^2=\\frac{25}{10}=2.5$\n",
    "2. $\\tau^2=2.56$​，因此$\\eta_n^2=\\frac{2.5 \\times 2.56}{2.5+2.56}=1.265$\n",
    "3. $\\mu=172.72, \\overline{x}=175.34$，因此$\\mu_n(\\overline{x})=\\frac{2.5 \\times 172.72+2.56 \\times 175.34}{2.5+2.56}=\\frac{431.8+448.8704}{5.06}=174.05$  \n",
    "所以$\\theta$的后验分布为：\n",
    "$$\n",
    "\\pi(\\theta|\\overline{x})=N(174.05,1.265) \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd4e207f",
   "metadata": {},
   "source": [
    "## 3 通过计算证明伽马分布和泊松分布是共轭分布\n",
    "\n",
    "结合韦书定义\"2.5.1\"：\n",
    "\n",
    "设$\\mathscr{F}$表示由$\\theta$的先验分布$\\pi(\\theta)$构成的分布族，如果对任取的$\\pi \\in \\mathscr{F}$及样本值$x$，后验分布$\\pi(\\theta|x)$仍属于$\\mathscr{F}$，则称$\\mathscr{F}$是一个共轭先验分布族。\n",
    "\n",
    "简单一点的理解共轭先验就是后验分布与先验分布属于同一分布族，只是分布函数中参数有变化。\n",
    "\n",
    "结合后验分布的计算公式：\n",
    "$$\n",
    "\\pi(\\theta|x)=\\frac{f(x|\\theta)\\pi(\\theta)}{\\int_\\Theta f(x|\\theta)\\pi(\\theta)\\mathrm{d}\\theta} \\tag{1}\n",
    "$$\n",
    "我们是可以用似然函数（likelihood function）来替代概率函数，即$l(\\theta|x)=f(x|\\theta)$。即：\n",
    "$$\n",
    "\\pi(\\theta|x)=\\frac{l(x|\\theta)\\pi(\\theta)}{\\int_\\Theta l(x|\\theta)\\pi(\\theta)\\mathrm{d}\\theta} \\tag{2}\n",
    "$$\n",
    "如果先验分布和似然函数可以使得先验分布和后验分布有相同的形式（属于同一分布族），那么就称先验分布与似然函数**是共轭的。**所以，共轭是指的先验分布和似然函数。\n",
    "\n",
    "下面来证明一下，伽马分布和泊松分布是共轭分布：\n",
    "\n",
    "我们假设总体服从泊松分布，$X=(X_1,X_2,...,X_n)$为从泊松分布$P(\\theta)$中抽取的随机样本，则样本$X$的概率函数为：\n",
    "$$\n",
    "\\begin{align*}\n",
    "f(x|\\theta)=l(x|\\theta)&=P(X_1=x_1,...,X_n=x_n|\\theta) \\\\\n",
    "&=\\prod_{i=1}^n\\frac{\\theta^{x_i}e^{-\\theta}}{x_i!}=\\frac{\\theta^{n\\overline{x}}e^{-n\\theta}}{x_1!...x_n!}\\propto \\theta^{n\\overline{x}}e^{-n\\theta} \\tag{3}\n",
    "\\end{align*}\n",
    "$$\n",
    "其中$\\overline{x}=\\frac{1}{n}\\sum_{i=1}^n x_i$表示样本的期望。\n",
    "\n",
    "我们令$\\theta$的先验分布是伽玛分布$\\Gamma(r,\\lambda)$，其密度函数为：\n",
    "$$\n",
    "\\pi(\\theta)=\\frac{\\lambda^r}{\\Gamma(r)}\\theta^{r-1}e^{-\\lambda\\theta}\\propto \\theta^{r-1}e^{-\\lambda\\theta} \\tag{4}\n",
    "$$\n",
    "将(4)式代入后验分布公式即公式(2)：\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\pi(\\theta|x)&\\propto l(x|\\theta)\\theta^{r-1}e^{-\\lambda\\theta} \\\\\n",
    "&\\propto \\theta^{n\\overline{x}}e^{-n\\theta}\\theta^{r-1}e^{-\\lambda\\theta} \\\\\n",
    "&\\propto \\theta^{n\\overline{x}+r-1}e^{-(n+\\lambda)\\theta}  \\tag{5}\n",
    "\\end{align*}\n",
    "$$\n",
    "将上式添加正则化因子后得到：\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\pi(\\theta|x)&=\\frac{(n+\\lambda)^{n\\overline{x}+r}}{\\Gamma(n\\overline{x}+r)}\\theta^{(n\\overline{x}+r)-1}e^{-(n+\\lambda)\\theta} \\\\\n",
    "&=\\Gamma(n\\overline{x}+r,n+\\lambda) \\tag{6}\n",
    "\\end{align*}\n",
    "$$\n",
    "这就说明，后验分布与先验分布属于同一分布族。即证明**伽马分布和泊松分布是共轭分布**。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70cd8721",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
